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We consider several models with long-range interactions evolving via Hamiltonian dynamics. The 
microcanonical dynamics of the basic Hamiltonian Mean Field (HMF) model and perturbed HMF 
models with either global anisotropy or an on-site potential are studied both analytically and nu- 
merically. We find that in the magnetic phase, the initial zero magnetization state remains stable 
above a critical energy and is unstable below it. In the dynamically stable state, these models 
exhibit relaxation time scales that increase algebraically with the number N of particles, indicat- 
ing the robustness of the quasistationary state seen in previous studies. In the unstable state, the 
corresponding time scale increases logarithmically in N. 

PACS numbers: 05.20.Gg,05.50.+q, 05.70.Fh 

I. INTRODUCTION 

In recent years, much work has been devoted to the understanding of the statistical mechanics and the dynamics of 
systems with long range interactions. In these systems the interaction potential between two particles decays at large 
distances r > 1 as V(r) ~ r~ a with a < d, d being the dimension of the system. Examples of such systems include 
sclf-gravitating systems, two-dimensional and geophysical vortices, non-neutral plasma and systems describing wave 
particle interactions (Free Electron Laser, CARL experiment etc.), and magnetic dipolar systems (see [l| for reviews). 

The long range nature of the interactions makes these systems non-additive. Due to this property, at statistical 
equilibrium, inequivalence between the microcanonical and the canonical ensembles is a generic feature. This was 
first observed in models of self gravitating stars [2, l3j] , and then seen in a number of other models ranging from the 
point vortex model 3, [p , plasma physics @, self gravitating systems @, Q> two dimensional flows @, long range 
Hamiltonian models [U[ll| to simple spin models with mean field interactions [H, [Hj]. A classification of phase 
transitions and ensemble inequivalence in generic long range systems [14| has shown that many possible types of 
behavior remain to be seen in specific physical systems. 

Beside these equilibrium peculiarities, the dynamics of systems with long range interactions also present several 
new features. For a large number N of particles, these systems may exhibit quasi-stationary states (QSS) [H, EH ( m 
the plasma or astrophysical context, see for instance UtI. [l8i) . very long relaxation time [10], vanishing Lyapounov 
exponents [l5|, [l|| , anomalous relaxation and diffusion [20l |21|. I22I [23l . l24j and breaking of ergodicity [TH, [25J . These 
features are a result of similar collective (self-consistent) dynamics [26| shared by systems with long range interactions. 
In the limit of large number of particles, such dynamics is well approximated by kinetic theories [H, HI, H3, HI, Hi! 
which to leading order in 1/y/N describe Vlasov type dynamics, and after much longer time, the relaxation towards 
equilibrium is governed by Lennard-Balescu type dynamics. 

In this paper, we consider the dynamics of systems with long range interactions and analyze the time it takes for a 
system to relax to its equilibrium state, starting from a thermodynamically unstable state. The relaxation process is 
usually initiated by the formation of droplets of the equilibrium state which is followed by a coarsening process [3(| • 
In systems with short range interactions, the initial droplets are of a typical radius which does not grow with the 
system size. As a result, the characteristic time for the formation of such droplets is also independent of the system 
size. On the other hand, in systems with long range interactions the relaxation time scale diverges with the system 
size N. This may result in long lived, quasistationary states (QSS) which are not the true equilibrium state of the 
system but which relaxes to the equilibrium state on time scales that increase algebraically with N. 

Such long-lived states have been seen, for instance, in numerical studies of long-ranged spin model [ill ]. These 
QS S have been explained as stable stationary states of the Vlasov equation and may lead to anomalous diffusion 
[23l . l3l| . We note that an alternative explanation, both for the existence of QSS and for anomalous diffusion has been 
proposed in the context of Tsallis non extensive statistical mechanics [H, |33[ (see [Hj] and (23[ for further discussions) . 
Recent studies have considered the possible prediction of QSS using the equilibrium statistical mechanics of the Vlasov 
equation [34], [H, H(| . The issue of the robustness of QSS when the Hamiltonian is perturbed by short range interactions 
[37l ] or when the system is coupled to an external bath [H, [39| has also been addressed, and it was found that while 
the power law behavior survives, the exponent may not be universal. Such a slow relaxation is not the only possible 
behavior in these systems [40| . A recent consideration of the thermodynamic stability of a mean field Ising model 
with stochastic dynamics has found the relaxation time to be logarithmic in N 13]. It is thus of interest to study 
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the slow relaxation processes in systems with long range interactions and explore in more detail the possible resulting 
time scales involved. 

With this aim in mind, we consider the microcanonical dynamics of a generalized Hamiltonian Mean Field (HMF) 
model which is a simple prototype of long-ranged systems by adding a global anisotro py te rm or an on-site potential 
energy term to the Hamiltonian. The basic HMF model ([4lJ, see also [HI, 0, Eallsj]) describes a system of N 
classical XY rotors with mean field coupling. Adding new terms to the Hamiltonian allows one to vary an external 
parameter (such as anisotropy) and explore a richer phase diagram. Our analysis of the Vlasov equation for these 
models shows that both logarithmic and power law behavior are generically present in long-ranged systems. It is found 
that at low energies the non-magnetic solution is dynamically unstable and the system relaxes to the magnetically 
ordered state on a logarithmic time scale which follows from the dynamic instability of Vlasov equation. At higher 
energies, but still within the magnetic state, the non-magnetic solution becomes linearly stable (although it is not the 
true equilibrium state of the system) and the relaxation takes place on algebraically diverging time scales and QSS 
are observed. We show the existence of QSS using analytic relations for the marginal stability of the Vlasov equation. 
These results give further insight into the robustness of QSS states when the interaction potential is perturbed. 

Most of the QSS studied so far have dealt with homogeneous situations (namely states whose distribution functions 
do not depend on the angle or spatial variable) . From a theoretical point of view, inhomogeneous QSS should exist in 
the same way as homogeneous ones. The main reason why such states have not been studied in detail is the difficulty 
to deal theoretically with the inhomogeneous marginal stability equation. The QSS of the HMF model with on-site 
potential described in Section llVI is an example of inhomogeneous QSS. Moreover, the new method we propose to 
study them can be applied to other situations, and also for instance to the usual isotropic HMF model. 

The relaxation times mentioned above are interesting from another viewpoint. Since the rationale for the existence 
of QSS is based on the approximation of the A particle dynamics by a Vlasov dynamics, a crucial issue is the 
understanding of the validity of such an approximation. In a classical work by Braun and Hepp [47| , it was proved 
that this approximation is valid for time t smaller than ty(N) for smooth interaction potential and large enough N. 
Following the reasoning of [43] allows one to conclude that a lower bound for ty (N) is proportional to In A for large 
A. A recent paper [48[ showed that, for initial conditions close to some homogeneous QSS, a lower bound for ty (A) 
scales like A 1 / 8 . In Section [V] we give two new results concerning this issue. First, we prove that the Braun and 
Hepp result is actually optimal. More precisely, we show that for some initial conditions, the kinetic description is 
not valid for time larger than ~ In N. Second, we argue that most of the trajectories will have ty either equal to 
the life time of a QSS (possibly algebraic), or logarithmically long ty depending on the way the relaxation towards 
equilibrium takes place. 

The rest of the paper is organised as follows. In Section [TH we define the HMF model and study its dynamically 
unstable state in detail using Vlasov equation and numerical simulations. The anisotropic HMF model is the subject 
of Section IIIII in which the dynamical phase diagram is obtained analytically in the energy-anisotropy plane. The 
model with on-site potential is introduced and discussed in Section ITVl The issue of the time scale over which Vlasov 
equation holds is discussed in Section [Vj Finally, we conclude with a summary and open questions in Section IVT1 



In this section, we consider the dynamics of the Hamiltonian Mean Field (HMF) model which is defined by the 
Hamiltonian 



where 9i and pi are the phase and momentum of the ith particle respectively, and A is the number of particles. 
In the equilibrium state, a second order phase transition between the ferromagnetic and paramagnetic state occurs 
at the critical energy density e c = 3/4. This has been shown in the canonical ensemble [4l[ and later verified for 
microcanonical ensemble using large deviations method (49j . 

The time evolution of this system starting far from the equilibrium state has been studied using Hamiltonian 
dynamics. The angle 9i and momentum pi of the z'th particle obey 



II. ISOTROPIC HAMILTONIAN MEAN FIELD MODEL 




(1) 



dt 
dpi 
dt 



—m x sin 9i + m y cos 9 t 



(2) 



(3) 
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where m x and m y are the components of the magnetization density 

(4) 



(^E cos ^E sin ^) ■ 

\ i=l i=l / 



The dynamics conserves the total energy and the total momentum. We start with an initial condition with randomly 
distributed 9i G [—tt, tt] so that the average magnetization is zero and the standard deviation about the mean is of 
the order ~ 1/y/N. To fix the total energy density e, the momentum is chosen from a distribution f^°'{p) with pi 
lying in the interval [p m i n , p m ax]- In this article, we consider the following choices of momentum distribution 



/ (l/2po) i P € [-po,Po] 



1 [P> ~ \ JWF^ exp (-(3p 2 /2) lP g(-oo,oo) K) 

where the parameter po in the uniform (or waterbag) distribution and /? in the Gaussian case are related to the energy 
as 

Po = V6e - 3 , 

= l/(2e-l) . (6) 

The equations of motion @ and §3§ are integrated using symplectic fourth order integrator with time step dt = 0.1. 
The reference coordinate axes in which 9i is measured is specified by the initial magnetization. 

To study the dynamical behavior of magnetization, we first recall the classical computation of the Vlasov equation 
US ES HE HH, [Ell (see also [U for the HMF model). The probability density f d (9,p,t) which counts the number of 
particles with angle 9 and momentum p at time t can be written as 



1 N 

f d (9,p,t) = -J2m(t)-e)6( Pi (t)-p) . (7) 
1=1 

Taking the time derivative of both sides of the above equation and using the canonical equations of motion, one 
obtains 

df d , df d dVdf d 

-dr +p ^9-^9^p- = °> (8) 

where the average potential V(9,t) is given by 

V(9,t) = T maX dp' f d9'{\ - cos{9 - 9'))f d (9',p',t) . (9) 



Expanding ([7]) to leading orders in 1/yN (see below), we obtain the Vlasov equation obeyed by the (smooth) 
distribution f(9,p,t) for infinite N, 

df df dVdf 

It is easily verified that the initial condition with angles distributed uniformly and momentum chosen from an arbitrary 
(normalised) distribution f^°'(p) is in fact a stationary state of this equation. To deal with the finite N case, we treat 
the finiteness as a perturbation about the homogeneous stationary state of the infinite system, 

f(9,p,t) = ±-f( \p) + \f^(9,p,t) , (11) 
where, after linearization, the perturbed distribution /W is a solution of the integro-differential equation 

^^^-i'-W- J ¥d0'f^(0',p',t) S H0-0')=0. (12) 

Since the initial angles and momentum of the N particles are sampled according to the distribution /'°', the small 
parameter A is of order 1 / y/N. 
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We now study the linear dynamics about the distribution f^{p)/2n by considering the eigenmodes of (jl2|) . A first 
treatment of the linear stability of stationary solutions of the HMF equation can be found in [46[ (see also [26} and 
fli}). As is explained for instance in textbooks on plasma physics [27} . any arbitrary function f^(d,p,t) cannot be 
decomposed into the eigenmodes of the linearized equation (|12[) . However when unstable modes exist, after a short 
time the largest of them dominates the dynamics and therefore at sufficiently long times, the temporal behavior can 
be described by the eigenmodes. We define the Fourier modes fj. (p, w) (and the conjugate /i^) as the eigenmodes 
of HUD of type : 

fV(6,p,t) = f^\p^)e^ ke+ ^ . (13) 

Since the last term in the Vlasov equation (fT2")) involves only e ±l9 , the Fourier modes must have k = ±1. The 
coefficients f±± are then determined by 

/±(1-J±) = (15) 



Integrating over p on both sides, one gets 



where 



t ! r _ dp 9f (0) 



7± = L, dp /ii(p) ' J± = "2 jL ■ (16) 

The frequency w is thus found from the condition J± = 1. For initial distributions which are even for the variable 
p, this condition yields an equation in lu 2 . We now consider specific choices of momentum distribution f^°\p). For 
uniformly distributed initial momentum, the frequency determined using the condition J± = 1 works out to be [4l| 

For e > e* = 7/12, unstable modes do not exist and the Vlasov equation is linearly stable. It is however unstable for 
e < e* and the perturbation f^'(6,p,t) grows exponentially fast towards the equilibrium state. Setting lu 2 = — Q 2 for 
fi real, we have 

fW(9,p,t) = Af { 1 1) (9,p)e m (18) 

where A is a constant. As mentioned above, this time dependence is valid for times t >> I/O. To treat the behavior 
at short times, the finite N behavior in the initial condition must be taken into account [271 ]. 

The average magnetization along the x and y axes, which are the observables of interest, can be written as 

fir 

(m x {t),m y (t)) = I dp d6 (cos0,sin0) f(6,p,t) (19) 

J p m i n J 7V 

= A r'"" dp f d9 {cosO, sin9) / (1) (0,p,t) +0{\ 2 ) . (20) 



The magnitude of the average magnetization is given by m — m 2 + m 2 and grows as 

1 



e nt . (21) 
N 



The results of our simulations in Fig. [T^ show that after a transient, the magnetization grows exponentially, as 
expected on the basis of preceding equation. Since the order of magnitude of the constant A in (|T8]) is proportional to 
1/VN, the scaled magnetization y/Nm x in Fig. |T|b collapses into a single curve. The growth rate is also in agreement 
with SI obtained in (I17|) . The above perturbative analysis cannot hold at long times as the linearization of the Vlasov 
equation breaks down when the magnetization reaches a value of order one. A similar analysis can be carried out 
for Gaussian distributed initial momentum. In this case, we obtain that the modes exist only for e < e* = 3/4 
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FIG. 1: (a) Time evolution of the average magnetization m(i) in the unstable phase at e = 0.55 for two values of N. (b) Data 
collapse for the scaled magnetization y~Nm x (t) vs. t on the semi logarithmic scale. The slope Q = l/y/5 of the solid line is 
given by (|17[) . The data have been averaged over 200 histories for each N. 



and the eigenvalues of the modes are given by (37) with D = 0. Thus, (21) is obeyed in this case as well with the 
corresponding eigenfrequency. Using (|21|) , we see that the time scale on which the system acquires a finite m diverges 
as In TV. Although this is the same behavior as in the thermodynamically unstable phase of Ising model [l3j], the 
origin of the logarithmic time scale here is dynamical while it follows from an argument based on thermodynamics in 
the Ising case. 

So far we have discussed the dynamical instability. The above argument gives a IniV time scale for relaxation 
for e < e*. Statistical mechanics predicts that homogeneous states are thermodynamically unstable for all (allowed) 
values of e < e c = 3/4. For e* < e < e c , the distribution f^{p) is dynamically stable but thermodynamically unstable. 
The relaxation is then not due to a dynamical instability and one then observes a N 1 - 7 time scale (see [la]). 



III. ANISOTROPIC HAMILTONIAN MEAN FIELD MODEL 



In this section, we consider the dynamics of the anisotropic HMF model defined by the Hamiltonian 

-. 2 



N r? 1 N 



8 = 1 



9j)] ~ ^ 



N 



cos 0i 



(22) 



where the last term represents the energy due to a global anisotropy in the magnetization along the x-axis. At zero 
temperature, if D is positive, the equilibrium magnetization is along the x-axis. Similarly for D < 0, the internal 
energy is lowered when the magnetization is along the y direction. For simplicity we consider below the case D > 0. 
To study the equilibrium phase diagram, consider the partition function in the canonical ensemble, 



N 

Z = I Y[ dpidOie 



-0H _ 



Y[d6i 



exp 



(23) 



On using a Hubbard-Stratonovich transformation [4l| , the integrals over angles can be rewritten as 



N 



2tt/V(1 + D) 



db x db y exp 



2(3 2(3(1 + D) 



In 



dO e 



b x cos 0-\-b v sin 9 



(24) 
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FIG. 2: Dynamical phase diagram in the (e, D) plane for the initial momentum distribution chosen as Waterbag (a) and 
Gaussian (b). The solid line is the thermodynamic phase boundary D(e c ) and the broken line is the dynamical boundary 
D(e*). The dotted line gives the line where to' vanishes. The dashed dotted line is the first order line where the magnetization 
changes its direction. Note that in the Gaussian case, the dynamical and the thermodynamical boundaries coincide. 



where the double integral can be evaluated using the saddle point method for large N, leading to the free energy per 
particle 



■InZ 



-1 



N(i 2/3 



In 



2/3(1 + D) 2/3 



In 



d6 cos ^~*~k y s * n ^ 



(25) 



In the above expression, b x and b y are determined by maximizing Ze with respect to b x and b y , and are related to the 
equilibrium magnetization fn x and m y along the x and y axis respectively as b x = /3(1 + D)rn x and b y = (3fn y . As 
explained above, for D > 0, the system orders along the x-axis and the magnetization rn x is determined by 



HI;, 



d9 cos 9 e 

J— 7T 



(3(1+0)^ 



r d9e 

J — 7T 



&{x+Dym x cos 9 



(26) 



Close to the critical point, the above transcendental equation can be expanded in a Taylor series about zero magneti- 
zation and we obtain m 2 x = (8/3(1 +D) - 16)/((4 - /3(1 + D))(i 2 (l + D) 2 ) for D > 0. The inverse critical temperature 
at which magnetization vanishes is given by /3 C = 2/(1 + D). The critical energy e c can be calculated using the free 
energy expression above and we obtain 



dp 



3 + D 



(27) 



For the anisotropic Hamiltonian, the equations of motion are similar to those for the isotropic case, © and ([3]), 
except that m x is replaced by (1 + D)m x . The average potential appearing in the corresponding Vlasov equation (JTT 
is now given by 



V(6,t) 



dp' / d6»'[l-cos(6»-6i')-£ | cos6'cos6' / ]/( 6 ''^'^) 



(28) 



It can be checked that the homogeneous state (in 9) is a stationary state of the Vlasov equation for the anisotropic 
HMF model so that the distribution f(0,p,t) can be written as ([IT]) . The distribution f^> is now a solution of the 
following equation, 



9 f(i) i 9 f(o) 
'P- 



dp'd9'f {1 \d',p',t) [sin(0 - 6') + Dsh 



dt r 89 2tt dp 
As before, going to the Fourier space and picking the coefficient of e ±l6 , we obtain 







/+ = 



D 



D 



J± 



(29) 



(30) 
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FIG. 3: Temporal behavior of the magnetization in the anisotropic HMF model for uniformly distributed initial momentum, 
(a) The three curves in the unstable phase show scaled magnetization y/~Nm x (t) for e = 0.55 and D — 0.9. The slope of 
the solid lines is given by uj in (|33|l . (b) Evolution in the stable phase with e = 0.8 and D = 0.9 (top to bottom) for 
N = 500(100), 2000(50), 5000(50), 10000(5). The number of histories over which data are averaged is given in parenthesis. The 
scaled data is consistent with N 1 ' 7 scaling of the quasistationary life time. 



where I± and J± are given in (|16[) . For D ^ 0, the frequency uj is then determined through 

(1 + D)J+J- - (l + |) (J+ + J_) + 1 = (31) 

which is a bilinear equation unlike in the isotropic case. We now find the dynamical phase diagram and the frequency 
in the unstable phase for two choices of initial momentum distribution. 

Uniform distribution.- For uniformly distributed /'°'(p), the integrals J± appearing in (|31| can be readily done and 
we obtain a fourth order equation for the frequency u), 

4(uj 2 - pi) 2 + 2 (2 + L>) [uj 2 -pi) + (1 + D) = (32) 

with the following four solutions 

S- t (.-I + °) (33) 
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-' 2 = 6^e--j. (34) 

The frequency to 2 vanishes at e* = (7 + D)/12 while uj' 2 becomes zero at e'* = 7/12. Thus in the magnetically 
ordered phase, e < e c , three regions can be identified: a linearly stable region e* < e < e c where long relaxation time 
scales increasing as a power of N are expected, an unstable regime e'* < e < e* with a single mode of instability namely 
fl, and another unstable regime e < e'* with two modes of instability fl and fi' where ft 2 = —uj 2 and fl' 2 = —uj' 2 . 
Since Q > f2', the magnetization increases exponentially fast with rate ft. Thus, in the last two regimes, the relaxation 
time scales diverging as IniV are expected. The (e, D) phase diagram resulting from this analysis is given in Fig. [5^. 

Our numerical results for the time evolution of the magnetization in the unstable and the stable phases are shown 
in Fig. [31 For e < e* (Fig. [3^), data collapse of the curves for various system sizes is observed when the magnetization 
is scaled with a factor ^/N as in the last section. Thus, we again obtain In AT scaling for the relaxation time. The 
growth rate in the unstable phase is also in agreement with uj in (f3"3"]l . For e > e*(Fig. [5Jd), the magnetization stays 
close to its initial value ~ 1/vN for a long time which is consistent with N 1 - 7 scaling as for the basic HMF model. 
Gaussian distribution.- For initial momentum chosen from Gaussian distribution, the dynamics are always unstable 
and the frequency uj 2 = — il 2 , f2 real. Consider the integrals J± defined in (fTB]) . 
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p -/3p 2 /2 

1-Q?\I— I dp— 




(35) 
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FIG. 4: Semi-log plot of the magnetization as a function of time for Gaussian f (p) in the dynamically unstable, ferromagnetic 
phase for N — 5 x 10 , D — 4 and various e. The lines have a slope Q given by (|37|l , 

where we have used that the derivative of /^°^ is an odd function. The integral in the last term can be evaluated 
using the Schwinger trick, 

Since the dispersion relation (f3"Tj) is quadratic in J + it has two solutions namely J+ = 1/(1 + D) and J+ = 1 with 
respective frequencies f2 and Q' which obey 

The real frequencies f2, f2' vanish at e* = (3 + D)/4 and e'* = 3/4 respectively. The critical energy e*(> e'*) coincides 
with e c as there is no stable phase when the initial momentum is distributed according to a Gaussian distribution. 
One is thus left with two unstable regimes, one with a single unstable mode (e'* < e < e c ), and the other with two 
unstable modes (e < e'*). Since the left hand side of the above equations for £1 and Q' is a monotonically decreasing 
function lying between 1 and 0, one has f2 > f2' for D > 0. Thus we expect the growth rate of m x to be f2 for all 
energies below e c . The dynamical phase diagram corresponding to the case of Gaussian initial distribution is given in 
Fig. ^Bp. Our numerical results for the evolution of the magnetization verifying the conclusions presented above are 
shown in Fig. |4j 



(3{l + D) 



2 





(37) 
(38) 



IV. HAMILTONIAN MEAN FIELD MODEL WITH ON-SITE POTENTIAL 

We next consider the HMF model with cosine on-site potential whose Hamiltonian is 

N N N 

h = ^ J2p 2 i + w £ (1 " cos{di -°^ +w J2 cos2 ^ • ( 39 ) 

where the last term in the Hamiltonian gives the energy due to an on-site potential. For positive (negative) W, the 
steady state magnetization is along the y(a;)-axis and in the following, we assume W > 0. The equilibrium properties 
of this model can be calculated following the same procedure outlined in Section [TTJ The equilibrium magnetization 
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vfiy along the y axis is determined by the following equation, 

d6 Sinfl e @m,y sin0-/Wcos 2 9 
m y = p 7q e f3m v sin6-f3Wcos 2 9 ' 

J — 7T 

The magnetization vanishes at a critical temperature determined via f3~ l = {Io{ z ) + Ii( z ))/^h){ z ) where Iq and I± 
are the the modified Bessel functions of the first kind and the argument z — (3 c W/2. For positive W -C 1, the critical 
temperature at which m becomes zero is given as (2 + W)/4 to leading order in W . The critical energy e c is then 
given by 

1 1 W 3 5W 

^w c + 2 + jr~^-- (41) 

For the model defined by (j3"9"|) , the canonical equations of motion reads, 

dOi dpi 

~dt = = — ™k sinwi + rrij, cosw, + 2VFcosc'i sin^i . (42) 

As in the previous sections, for a large number of particles the dynamics is well approximated by the Vlasov equation 
(TTU)) with the following potential 

V[f] = - dp' [ d6'cos(6-d')f(6',p',t) + Wcos 2 6 . (43) 



Unlike for the models considered in the preceding sections, the homogeneous state is no longer a stationary state 
of the Vlasov equation due to the presence of the last term in the potential V. However, the distribution function 
f(6,p) = <f>(e(6,p)), with arbitrary function $, where the single particle energy e is given by 

2 2 

e (0,p) = ?^+V{9) = — — tti x cos 9 — my sin 9 + W cos 2 9 . (44) 

are stationary solutions of the Vlasov equation (fit))) . This can be easily seen by differentiating Vlasov equation (jTUJ) 
and by noting that the potential energy V is a function of the angle variable only. We stress that the magnetization 
values m x and m y must be self consistently determined. The particular case (9,p) ~ exp (— (3e (9,p)) is the 
statistical equilibrium density. 

A special class of stationary distributions is given by f^°'(9,p) = < &(p 2 /2 + VFcos 2 9). The fact that m x — m y — 
follows by symmetry arguments. In the following, we will study a waterbag stationary state distribution i.e. <!> is a 
step function 

A ,£ + Wcos 2 9<E 
, otherwise . 

The distribution function f^(p) is thus constant over the domain T> defined by 



ft°Xe, P ) = \ A , £ t + wcos^<e (45) 

, otherwise . 



\p\ < p Q (9) = y/2 (E-W cos 2 9) , (46) 

and is simply connected for W < E. The study of the case W > E can be done following ideas similar to those 
described below. The normalization constant A is determined using d9 Jfp^m dp f^(9,p) — 1 and we have 



- 1 - = / dd y / 2(E- Wcos 2 i 



(47) 



The parameter E can be related to the conserved initial energy e by performing the integration over the p variable in 
(I3"9")) and we obtain 



A J dd I - [2 (E - W cos 2 6»)] 3/2 + 2W cos 2 6 [2 (E - W cos 2 9)] 



1/2 



(48) 



To compute the linear stability threshold, as before, we linearize the dynamics close to the stationary solution, 
/ = / (0) + A/W exp(zwt). The perturbation f^- 1 ' then satisfies 

; +p d9 dp d9 dp ~ [ ' 
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We will determine the energy e* which corresponds to the neutral mode u> = 0. Below this energy, the dynamics is 
expected to be unstable and stable above it. At u> = 0, the above equation can be written explicitly as 

P-^q--A[8(p + pq(P))-6(p-po(9))] J dp'd6'sm{6-8')f( 1 \8',p',u) + 2Wsm0cos6^— = . (50) 

We solve this last equation by a formal expansion in terms of Dirac distributions 5 and its order n derivatives §^ n ' : 

oo 

= E a « W 5(n) (P + W) + 6 » W (P " W) ( 51 ) 

The equations for a„ and 6„ and their solutions can be found recursively. We report the analysis only for ao and bo : 

da Q _ db Q _ sin 9m x [Sf] - cos 9m y [Sf] 

d9 ~ d9 p„(0) l ° > 

The magnetization [5/] must be determined self-consistently from the distribution function. Using (|5ip . the 
formula (|52|) and (|47|. we obtain 

<*V2(£-WW0)= f dB-==^== . (53) 
J-tt (i? — W cos ^) 

This is the equation for the marginal stability of the inhomogeneous waterbag distribution function. We define 

m = 1 f *> (54) 

J 27r d9 ^l~xcos 2 9 Jo Vl-xcos 2 ^ 

which can be expressed in terms of complete elliptic functions of the first and second type. The equation for the 
marginal stability for inhomogeneous Water Bag distributions then reads 

2E = l(^j . (55) 

One can prove that I is a strictly decreasing function from the interval [0; 1] onto the interval [1/2; 1]. 7 _1 is thus 
an increasing function from [1/2; 1] onto [0; 1]. From this, one can prove that, the equation I^ 1 (2E) = W/E has a 
single solution E* (W), for each value of W in the range [0; 1/2], and no solutions for W > 1/2. For W — 1/2, we 
have E* = W ; this is the limit above which the inhomogeneous Water-Bag cease to be simply connected (see the 
discussion below (|1d]) ). For values larger than W = 1/2, the transition value E* (W) could be studied by considering 
doubly connected domains. 

For W « 1, equation (153)) can be easily linearized and we obtain E* (W) = 1/4 + 3W/8 + O (W 2 ). Now using 
we can compute the critical e* (W) as a function of W. For instance, for W << 1, we obtain 

e*(W) = ^ + ^ + 0(W 2 ) . (56) 

For W — 0, we obtain the critical energy for the homogeneous Water Bag distribution e* = 7/12 (see (3). For any 
value < W < 1/2, the critical value of the energy can easily be computed numerically from ([55]) . Figure [5] shows 
the curve of marginal stability. 

At energies below e*(W), where the waterbag state is unstable, the relaxation time grows logarithmically with the 
system size as in the case of the model with global anisotropy. This is a result of the fact that the initial magnetization 
of the waterbag state is of order For energies above e*(W) where the waterbag state is linearly stable, we 

expect a longer relaxation time, which grows algebraically with N. 

V. HOW LONG IS THE VLASOV APPROXIMATION VALID ? 

In this section, we discuss the validity of the approximation of the A^-particle phase-space distribution 

1 N 

fd(0,p, t) = - WW - e ) 5 (Pi(t) - P) (57) 
i=l 
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FIG. 5: Dynamical phase diagram in the (e, W) plane for the inhomogeneous uniform distribution of the HMF model with 
on-site potential. The bold line is the line of marginal stability W(e*) and the thermodynamic phase boundary W(e c ) is shown 
as a dashed line. The dotted lines are the respective curves obtained within linear approximation for W <C 1 given in (|4ip and 
(SB). 



by a smooth density function f s (6,p,t), where both fd and f s satisfy the Vlasov equation (fTO)) (mean field approxi- 
mation). We ask: how long is this approximation justified? 

We consider an ensemble of initial conditions for the iV-particles {(Oi(t = 0),p;(i = 0))i<i<jv}. We suppose that 
the corresponding fd{t = 0) is close to some smooth distribution function f s (t — 0) (i.e. the distance between fd 
and f s converges to zero as N goes to infinity). If for almost all initial conditions, fd remains close to the smooth 
f s , we say that the system has a kinetic behavior: all trajectories remain close to each other. The kinetic evolution 
is then the solution of the Vlasov equation with initial condition f s (t = 0). The issue of the validity of the mean 
field approximation, or equivalently of the validity of the kinetic description, is to know for how long all trajectories 
remain close to f s . Since fd satisfies the Vlasov equation, this issue is related to the stability of the Vlasov equation. 

Let us first summarize the known results. The validity of this mean field approximation for large N has been 
established mathematically, for smooth potential V, by Braun and Hepp [13] (see also [!(|). More precisely, the 
theorem of Braun and Hepp states that for a mean-field microscopic two-body smooth potential, the distance between 
two initially close solutions of the Vlasov equation increases at most exponentially in time. Indeed, if f(t) and g(t) are 
two solutions, if A(f,g)(t = 0) is sufficiently small, then A(f,g)(t) < A(/, <?)(0) exp(at), where A is the Wasserstein 
distance and a is a constant. 

This result can be applied to the approximation of the TV-particle Hamiltonian dynamics by the Vlasov equation. 
We consider, for the N particle dynamics, an ensemble of initial conditions {(#i,Pi)i<i<Ar} distributed according to 
the measure f(0±,p±, ...^OmtPn) — ^-i—xfsiPi^Pi)- F° r large N, for a typical initial condition, the phase-space initial 
distribution f d (0,p,0) will be close to f s (6,p,0). Typically A(f d ,f s )(t = 0) = 0(1/N a ) with a > 0. Let us consider 
f s the solution of the Vlasov equation with initial condition /W. Because f d and f s are both solutions for the Vlasov 
equation, we can apply the Braun and Hepp's result. If we define ty to be the time at which the error A(fd, f s ) 
is of order unity, the theorem then implies that ty increases at least as lnJV when N — + oo (for this argument, see 
also 

In Sections HT1 and HTT1 we have considered the special case where the initial distribution f^ is a stationary solution 
of the Vlasov equation. When /(°' is unstable, the perturbation /W grows exponentially, we have explained and 
illustrated that ty is proportional to In TV. Thus the trajectories diverge from f s on a time scale given by InTV/fL 
This thus proves that the Braun and Hepp's result for ty is not only a lower bound, but is actually achieved. After 
this time scale, the trajectories diverge and the system does not have a kinetic behavior any more. 

When /(°) is a stable stationary solution of the Vlasov equation, ty is the stability time of the Quasi-Stationary 
State. A very recent work [48| has proven that the A^-particle dynamics actually remains close to the stationary 
solutions, at least for times of order N 1 / 8 , when the potential V is sufficiently smooth. On physical grounds, using 
kinetic theory, one expects the validity time to be of order N for systems in which each particle is characterised by 
more than one dynamical variable. When no resonance between trajectories is possible as is the case for systems with 
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one dynamical variable, one expects validity times to be much larger than N [23j] . This peculiarity of Id systems 
has been numerically observed in the HMF model, where times of order iV 17 have been measured for homogeneous 
Quasi-Stationary states 16]. This time scale seems to be robust with respect to the perturbations of the Hamiltonian 
as we have demonstrated for the anisotropic case in Section IIIII 

Let us now consider the more general case when the initial condition is close to a distribution / which is not 
stationary. In accordance with the observed phenomenology for the Vlasov equation, one expects that / will have a 
rapid relaxation in a finite time either towards a Quasi-Stationary state, or towards a periodic solution, or towards 
a statistical equilibrium for the Vlasov equation. Because this first stage takes place in times which are of order one 
(which do not depend on N), one expects that this initial relaxation will have a negligible effect on the long time error. 
The validity time for the Vlasov approximation will then be given by the validity time of the Quasi-Stationary states. 
Then times ty of order N a are expected. We thus conjecture that, for generic initial distributions, the approximation 
by the Vlasov equation is valid over times which are the life time for the Quasi-Stationary states. 

VI. CONCLUSIONS 

In this article, we studied the short time dynamics of models with long-ranged Hamiltonians. In each case, starting 
from the initial magnetization of order 1 / \/N, the time required to achieve a finite value of m scales with the number 
N of particles. This behavior is different from that of the corresponding short-ranged Hamiltonians where such time 
is of order unity. The dynamics were studied by a numerical integration of the Hamilton's equations of motion and 
a stability analysis of Vlasov equation. The latter analysis shows that close to the unstable stationary states of the 
Valsov equation, the relaxation occurs over a time ~ In N while close to stable stationary states, the system stays in 
a quasistationary state whose lifetime goes as a power law in N. 

So far we considered only the deterministic dynamics. An interesting direction would be to study such models 
when the dynamical rules are stochastic. One such case has been discussed in [l3[ where the system of Ising spins 
evolves via microcanonical Monte Carlo dynamics. Since the Ising energy can be obtained as a limiting case of the 
Hamiltonians considered here (with anisotropy or on-site potential), it would be worthwhile to study these systems 
with stochastic evolution rules. Besides, for the HMF model, the ensemble equivalence in the equilibrium steady 
state has been shown using large deviation method [49[ . Recent numerical studies [H, [H| of the dynamics of this 
system in contact with a thermal bath have focused on the stable regime e > e* and find that the quasistationary 
states seen in the microcanonical ensemble survive but the lifetime increases as a power law (in N) with an exponent 
that decreases with increasing system-bath coupling. It would be interesting to know if the unstable phase exhibits a 
similar dependence on the coupling with the heat reservoir. 

Acknowledgement: Financial support of the Israel Science Foundation (ISF) is acknowledged. Visit of F.B. to the 
Weizmann Institute has been supported by the Albert Einstein Minerva Center for Theoretical Physics. 



[1] Dynamics and Thermodynamics of Systems with Long Range Interactions, T. Dauxois, S. Ruffo, E. Arimondo, M. Wilkens 

Eds., Lecture Notes in Physics Vol. 602, Springer (2002). 
[2] D. Lynden-Bell and R. Wood, Mon. Not. R. Astron. Soc.138 495 (1968). 
[3] P. Hertel and W. Thirring, Annals of Physics 63, 520 (1971). 
[4] R.A. Smith and T.M. O'Neil, Phys. Fluid. B 2, 2961 (1990). 

[5] E. Caglioti, PL. Lions, C. Marchioro and M. Pulvirenti, Commun. Math. Phys. 174 229 (1990). 
[6] M.K.H. Kiessling and T. Neukirch, Proc. Nat. Acad. Sci. USA 100, 1510 (2003). 
[7] B.N. Miller and P. Youngkins, Phys. Rev. Lett. 81, 4794 (1998). 
[8] P.H. Chavanis and I. Ispolatov, Phys. Rev. E 66, 036109 (2002). 
[9] R.S. Ellis, K. Haven, and B. Turkington, Nonlinearity 15, 239 (2002). 
[10] M. Antoni, S. Ruffo and A. Torcini, Phys. Rev. E 66, 025103 (2002). 

[11] T. Dauxois, S. Lepri and S. Ruffo, Commun. Nonlinear Sci. Numer. Simul. 8, 375 (2003). 

[12] J. Barre, D. Mukamel, S. Ruffo, Phys. Rev. Lett. 87, 030601 (2001). 

[13] D. Mukamel, S. Ruffo and N. Schreiber, Phys. Rev. Lett. 95, 240604 (2005). 

[14] F. Bouchet, J. Barre, J. Stat. Phys. 118, 1073 (2005). 

[15] V. Latora, A. Rapisarda, S. Ruffo, Phys. Rev. Lett. 80, 692 (1998). 

[16] Y. Y. Yamaguchi, J. Barre, F. Bouchet, T. Dauxois and S. Ruffo, Physica A 337, 36 (2004). 
[17] D.H.E. Dubin and T.M. O'Neil, Rev. Mod. Phys. 71, 1 (1999). 

[18] L. Spitzer, Dynamical evolution of globular clusters, Princeton Series in Astrophysics (1987). 
[19] M.C. Firpo, Phys. Rev. E 57, 6599 (1998). 

[20] V. Latora, A. Rapisarda, S. Ruffo, Phys. Rev. Lett. 83, 2104 (1999). 



13 



Y.Y. Yamaguchi, Phys. Rev. E 68, 066210 (2003). 

A. Pluchino, V. Latora, A. Rapisarda, Phys. Rev. E 69, 056113 (2004). 

F. Bouchet and T. Dauxois, Phys. Rev. E 72, 045103 (2005). 

Y. Yamaguchi, F. Bouchet and T. Dauxois, J. Stat. Mech. 1 P01020 (2007). 

F. Borgonovi, G. L. Celardo, M. Maianti, E. Pedersoli, J. Stat. Phys. 116, 1435 (2004). 

D. Del-Castillo-Negrete, Physica A 280, 10 (2000). 

D.R. Nicholson, Introduction to plasma theory, Wiley New York (1983). 
D.H.E. Dubin, Physics of Plasmas 10, 1338 (2003). 
P.H. Chavanis, Phys. Rev. Lett. 84, 5512 (2000). 
A.J. Bray, Adv. Phys. 43, 357 (1994). 

F. Bouchet, T. Dauxois, Journal of Physics: Conference Series 7, 34 (2005). 
V. Latora, A. Rapisarda, C. Tsallis, Phys. Rev. E 64, 056134 (2001). 
A. Pluchino, V. Latora, A. Rapisarda, Physica A 340, 187 (2004). 

J. Barre, T. Dauxois, G. De Ninno, D. Fanelli and S. Ruffo, Phys. Rev. E 69, 045501(R) (2004). 

A. Antoniazzi, D. Fanelli, J. Barre, T. Dauxois. P.H. Chavanis and S. Ruffo, Phys. Rev. E 75, 011112 (2006). 

P.H. Chavanis, Eur. Phys. J. B 53, 487 (2006). 

A. Campa, A. Giansanti, D. Mukamel and S. Ruffo, Physica A 365, 120 (2006). 
F. Baldovin, E. Orlandini, Phys. Rev. Lett. 96, 240602 (2006). 

F. Baldovin, E. Orlandini, Phys. Rev. Lett. 97, 100601 (2006). 

Relaxation times can be exponentially long if the initial state is metastable which we do not consider here. 
M. Antoni and S. Ruffo, Phys. Rev. E 52, 2361 (1995). 

G. M. Zaslavsky, V. F. Shabanov, K. S. Aleksandrov, I. P. Aleksandrov, Sov. Phys. JETP 45, 315 (1977). 
T. Konishi and K. Kaneko, J. Phys. A 25, 6283 (1992). 

D. del-Castillo-Negrete, in Dynamics and Thermodynamics of Systems with Long Range Interactions, T. Dauxois, S. Ruffo, 

E. Arimondo, M. Wilkens Eds., Lecture Notes in Physics Vol. 602, Springer (2002). 
C. Pichon, PhD thesis, Cambridge (1994). 

S. Inagaki and T. Konishi, Publ. Astron. Soc. Japan 45, 733 (1993). 

W. Braun and K. Hepp, Comm. Math. Phys. 56, 101 (1977). 

E. Caglioti and F. Rousset, Comm. Math. Sci. 5, 11 (2007). 

J. Barre, F. Bouchet, T. Dauxois and S. Ruffo, J. Stat. Phys. 119, 677 (2005). 

H. Spohn, Large scale dynamics of interacting particles, Springer- Verlag New York (1991). 
M.C. Firpo and Y. Elskens, J. Stat. Phys. 93, 193 (1998). 

Y. Elskens and D.F. Escande, Microscopic dynamics of plasmas and chaos, IoP Publishing, Bristol (2002). 
H. Morita and K. Kaneko, Europhys. Lett. 66, 198 (2004). 



